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ABSTRACT 

The intense turbulence present in the solar convection zone is a major challenge to both theory 
and simulation as one tries to understand the origins of the striking differential rotation profile 
with radius and latitude that has been revealed by helioseismology. The differential rotation 
must be an essential element in the operation of the solar magnetic dynamo and its cycles of 
activity, yet there are many aspects of the interplay between convection, rotation and magnetic 
fields that are still unclear. We have here carried out a series of 3-D numerical simulations of 
turbulent convection within deep spherical shells using our anelastic spherical harmonic (ASH) 
code on massively parallel supercomputers. These studies of the global dynamics of the solar 
convection zone concentrate on how the differential rotation and meridional circulation are 
established. We have addressed two issues raised by previous simulations with ASH. Firstly, 
can solutions be obtained which possess the apparent solar property that the angular velocity 
f2 continues to decrease significantly with latitude as the pole is approached? Prior simulations 
had most of their rotational slowing with latitude confined to the interval from the equator to 
about 45°. Secondly, can a strong latitudinal angular velocity contrast Afi be sustained as the 
convection becomes increasingly more complex and turbulent? There was a tendency for Afi 
to diminish in some of the turbulent solutions that also required the emerging energy flux to 
be invariant with latitude. 

In responding to these questions, five cases of increasingly turbulent convection coupled with 
rotation have been studied along two paths in parameter space. We have achieved in one case 
the slow pole behavior comparable to that deduced from helioseismology, and have retained 
in our more turbulent simulations a consistently strong Afi. We have analyzed the transport 
of angular momentum in establishing such differential rotation, and clarified the roles played 
by Reynolds stresses and the meridional circulation in this process. We have found that the 
Reynolds stresses are crucial in transporting angular momentum toward the equator. The effects 
of baroclinicity (thermal wind) have been found to have a modest role in the resulting mean 
zonal flows. The simulations have produced differential rotation profiles within the bulk of the 
convection zone that make reasonable contact with ones inferred from helioseismic inversions, 
namely possessing a fast equator, an angular velocity difference of about 30% from equator to 
pole, and some constancy along radial lines at mid-latitudes. Future studies must address the 
implications of the tachocline at the base of the convection zone, and the near-surface shear 
layer, upon that differential rotation. 

Subject headings: convection - hydrodynamics - Sun: interior - Sun: rotation - turbulence 



1. INTRODUCTION 

The solar turbulent convection zone has striking 
dynamical properties that continue to challenge ba- 
sic theory. The most fundamental issues involve the 
solar rotation profile with latitude and depth, and the 
manner in which the 22-year cycles of solar magnetic 
activity are achieved. These two issues are closely in- 



terrelated, for the global dynamo action is likely to be 
very sensitive to the angular velocity Vl profiles real- 
ized by convection redistributing angular momentum 
within the deep zone. Both dynamical topics touch 
on the seeming inconsistency that turbulence can be 
both highly intermittent and chaotic on smaller spa- 
tial and temporal scales, yet exhibit large-scale or- 
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dered behavior (cf. Brummell, Cattaneo & Toomre 
1995). The differential rotation profile established 
by the turbulent convection, though strong in con- 
trast, is remarkably smooth; the global-scale mag- 
netic activity is orderly, involving sunspot eruptions 
with very well defined rules for field parity and emer- 
gence latitudes as the cycle evolves. The wide range 
of dynamical scales of turbulence present in the solar 
convection zone yield severe challenges to both the- 
ory and simulation: the discernible structures range 
from granules (~ 10 3 km or 1 Mm in horizontal size), 
to supergranules (~30 Mm), to possible patterns of 
giant cells comparable to the overall depth of that 
zone (~200 Mm, or nearly 30% by radius). Given 
that the dissipation scales are on the order of 0.1 
km or smaller, the solar turbulence encompasses at 
least six orders of magnitude for each of the three 
physical dimensions. The largest current 3-D turbu- 
lence simulations can resolve about three orders of 
magnitude in each dimension. Yet despite the vast 
difference in the range of scales dynamically active in 
the sun and those accessible to simulations, the lat- 
ter have begun to reveal basic self-ordering dynamical 
processes yielding coherent structures that appear to 
play a crucial role in the global differential rotation 
and magnetic dynamo activity realized in the sun. 

It has long been known by tracking surface features 
that the surface of the sun rotates differentially (e.g. 
Ward 1966, Schiissler 1987): there is a smooth pole- 
ward decline in the angular velocity f2, the rotation 
period being about 25 days in equatorial regions and 
about 33 days near the poles. Helioseismology, which 
involves the study of the acoustic p-mode oscillations 
of the solar interior (e.g. Gough & Toomre 1991), 
has provided a remarkable new window for studying 
dynamical processes deep within the sun. This has 
been enabled by nearly continuous helioseismic ob- 
servations provided from the SOHO spacecraft with 
the high-resolution Michelson Doppler Imager (SOI- 
MDI) (Scherrer et al. 1995) and from the ground- 
based Global Oscillation Network Group (GONG) 
set of six related instruments (Harvey et al. 1996). 
The helioseismic findings about differential rotation 
deeper within the sun have turned out to be revolu- 
tionary, for they are unlike any anticipated by con- 
vection theory prior to such probing of the interior of 
a star. Helioseismology has revealed that the rotation 
profiles obtained by inversion of frequency splittings 
of the p modes (e.g. Libbrecht 1989, Thompson et al. 
1996, Schou et al. 1998, Howe et al. 2000b) have the 
striking behavior shown in Figure |l[ The variation of 
angular velocity f2 observed near the surface, where 
the rotation is considerably faster at the equator than 
near the poles, extends through much of the con- 
vection zone with relatively little radial dependence. 



Thus at mid-latitudes f2 is nearly constant on radial 
lines, in sharp contrast to early numerical simulations 
of rotating convection in spherical shells (e.g. Gilman 
& Miller 1981, Glatzmaier 1987) that suggested that 
f2 should be nearly constant on cylinders aligned with 
the rotation axis and decreasing inward on the equa- 
torial plane. Another striking feature is the region of 
strong shear at the base of the convection zone, now 
known as the tachocline, where adjusts to apparent 
solid body rotation in the deeper radiative interior. 
Whereas the convection zone exhibits prominent dif- 
ferential rotation, the deeper radiative interior does 
not; these two regions are joined by the complex shear 
of the tachocline. There is further a thin shear bound- 
ary layer near the surface in which £1 increases with 
depth at intermediate and high latitudes. 

The tachocline has been one of the most surpris- 
ing discoveries of helioseismology, especially since its 
strong rotational shear affords a promising site for the 
solar global dynamo. Such a tachocline was not antic- 
ipated, and current theoretical approaches to explain 
its presence are still only innovative sketches (Spiegel 
& Zahn 1992; Gough & Mclntyre 1998; Charbon- 
neau, Dikpati & Gilman 1999). Helioseismology has 
also recently detected prominent variations in the ro- 
tation rate near the base of the convective envelope, 
with a period of 1.3 years evident at low latitudes 
(Howe et al. 2000a; Toomre et al. 2000). These 
are the first indications of dynamical changes close 
to the presumed site of the global dynamo as the 
cycle advances. Such a succession of developments 
from helioseismology provide both a challenge and a 
stimulus to theoretical work on solar convection zone 
dynamics. 

Seeking to understand solar differential rotation 
and magnetism requires 3-D simulations of convec- 
tion in the correct full spherical geometry. However, 
the global nature of such solutions represent a major 
computational problem, given that the largest scale 
is pinned and only three orders of magnitude smaller 
in scale can be represented. Much of the small-scale 
dynamics in the sun dealing with supergranulation 
and granulation are by necessity then largely omit- 
ted. The alternative is to reduce the fixed maximum 
scale by studying smaller localized domains within 
the full shell and utilizing the three orders of mag- 
nitude to encompass the dynamical range of turbu- 
lent scales. There are clear tradeoffs: the global 
models operate in the correct geometry yet struggle 
to encompass enough of a dynamical range to ad- 
mit fully turbulent solutions, whereas the local mod- 
els are able to study intensely turbulent convection 
but only within a particular limited portion of the 
full domain. Both approaches are needed, and the 
efforts are complementary, as reviewed in detail by 
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Fig. 1. — (a) Angular velocity profile Q./2n with radius and latitude as deduced from helioseismology using SOI-MDI data, 
with red tones indicating fast rotation and blue-green the slowest rotation [adapted from Schou et al. 1998]. (b) Time-averaged 
rotation rates from five years of GONG helioseismic data, plotted against radius at different latitudes. The surface shear layer and 
the tachocline at the base of the convective zone are indicated, as well as the zone covered by our computational domain (grey 
area) [adapted from Howe et al. 2000b]. 



Gilman (2000) and Miesch (2000). Highly turbulent 
but localized 3-D portions of a convecting spherical 
shells are being studied to assess transport properties 
and topologies of dynamical structures (e.g. Bran- 
denburg et al. 1996; Brummell et al. 1996, 1998; 
Porter & Woodward 2000; Robinson & Chan 2001), 
of penetration into stable domains below (Brummell 
et al. 2001, Porter & Woodward 2001), of effects of 
realistic near-surface physics on granulation and su- 
pergranulation (e.g. Stein & Nordlund 1998), and 
of dynamo processes and magnetic transport by the 
convection (e.g. Cattaneo 1999, Tobias et al. 2001). 
Without recourse to direct simulations, the angular 
momentum and energy transport properties of tur- 
bulent convection have also been considered using 
mean-field approaches to derive second-order corre- 
lations (the Reynolds stresses and anisotropic heat 
transport) under the assumption of separability of 
scales. Although such procedures involve major un- 
certainties, the resulting angular momentum trans- 
port, which is described by mechanisms such as the 
so-called A effect, have served to reproduce the so- 
lar meridional circulation (e.g. Durney 1999, 2000) 
and differential rotation (e.g. Kichatinov & Riidiger 
1995). Various other states can be achieved by ad- 
justing parameters. 

Initial studies of convection in full spherical shells 
to assess effects of rotation with correct account of 
geometry (e.g. Gilman & Miller 1981; Glatzmaier & 
Gilman 1982; Glatzmaier 1985, 1987; Sun & Schubert 
1995) have set the stage for our efforts to study more 
turbulent flows using new numerical codes designed 



for the massively parallel computer architectures that 
are enabling such major simulations. We here report 
on our continuing studies with the anelastic spherical 
harmonic (ASH) code (Clune et al. 1999) to examine 
the O, profiles established within the bulk of the solar 
convection zone by turbulent convection, building on 
the progenitor work by Miesch et al. (2000), Elliott, 
Miesch & Toomre (2000), and Brun & Toomre (2001). 
We also recognize the recent modelling of convection 
in spherical shells by Takehiro & Hayashi (1999) and 
Grote & Busse (2001). 

The simulations reported in Miesch et al. (2000) 
and Elliott et al. (2000) have revealed the richness 
and complexity of compressible convection achieved 
in rotating spherical shells. Most of the resulting an- 
gular velocity profiles in the seven simulations con- 
sidered have begun to make substantial contact with 
the helioseismic deductions within the bulk of the so- 
lar convection zone. These possess fast equatorial 
rotation (prograde), substantial O contrasts with lat- 
itude, and reduced tendencies for rotation to be con- 
stant on cylinders. The simulations with ASH have 
not yet sought to deal with questions of the near- 
surface rotational shear layer nor with the forma- 
tion of a tachocline near the base of the convection 
zone. These studies have revealed that to achieve 
fast equators it is essential that parameter ranges be 
considered in which the convection senses strongly 
the effects of rotation, which translates into having 
a convective Rossby number less than unity for large 
Taylor numbers. Such rotationally-constrained con- 
vection exhibits downflowing plumes that are tilted 
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away from the local radial direction, resulting in ve- 
locity correlations and thus Reynolds stresses that are 
found to have a significant role in the redistribution 
of angular momentum. This seems to provide paths 
to realize solar-like ft profiles. Further, it is desirable 
to impose thermal boundary conditions at the top of 
the domain that enforce the constancy of emerging 
flux with latitude in order to be consistent with what 
appears to be observed. 

We wish to focus on two outstanding issues raised 
by the prior simulations with ASH that need par- 
ticular attention concerning the differential rotation 
established within the bulk of the solar convection 
zone. As Issue 1, the helioseismic inferences in Fig- 
ure 1 emphasize that ft in the sun appears to decrease 
significantly with latitude even at mid and high lati- 
tudes, a property which has been difficult to attain in 
the prior seven simulations. The substantial latitudi- 
nal decrease in angular velocity, say Aft, in the mod- 
els is primarily achieved in going from the equator to 
about 45°, with little further decrease in ft achieved 
at higher latitudes in most of the cases. Whereas the 
overall latitudinal contrasts from equator to pole in 
the models and in the sun are roughly of the same or- 
der, the angular velocity in the sun continues to slow 
down much more as the pole is approached. Two 
models, designated as LAM (in Miesch et al. 2000) 
and L3 (in Elliott et al. 2000), do exhibit ft which de- 
crease at high latitudes, but LAM involves an emerg- 
ing heat flux that varies too much with latitude due to 
choice of boundary conditions, and L3 has an overall 
Aft that is only two-thirds of the helioseismic value. 
Thus in confronting Issue 1, we will search in param- 
eter space for solutions that can achieve ft profiles in 
which the decrease with latitude does not taper off 
at mid latitudes and for which the contrast Aft is at 
least comparable to the helioseismic findings. 

As Issue 2, with the convection becoming more tur- 
bulent, achieved by decreasing either the thermal or 
viscous diffusivities, there is a tendency for the latitu- 
dinal contrast Aft in the solutions to diminish or even 
decrease very prominently, thus being at variance 
with Aft deduced from helioseismology. This behav- 
ior appears to arise from increasing complexity lead- 
ing to a weakening of nonlinear velocity correlations 
that have a crucial role in angular momentum redis- 
tribution. These Reynolds stress terms are strong 
in the laminar solutions that involve tilted colum- 
nar convection cells ('banana cells') aligned with the 
rotation axis; they weaken as the flows become more 
intricate, but would be expected to become again sig- 
nificant once coherent structures develop at higher 
levels of turbulence. For example, the model TUR 
(in Miesch et al. 2000) exhibits the emergence of 
downflow networks involving fairly persistent plumes 



which possess some of the expected attributes of the 
coherent structure seen in localized domains of highly 
turbulent convection (e.g. Brummell et al. 1998). As 
a result, TUR has a fairly interesting angular mo- 
mentum transport attributed to the nonlinear cor- 
relations that sustain a level of differential rotation 
slightly weaker than LAM, but it too has a consider- 
able variation of heat flux with latitude. The model 
T2 (in Elliott et al. 2000) sought to correct the lat- 
ter by using modified thermal boundary conditions 
but appears to not have attained high enough tur- 
bulence levels to realize strong coherent structures. 
Absent those features, T2 yielded ft profiles with a 
small Aft, and even a slightly slower equatorial ro- 
tation rate than that in the mid latitudes. Thus 
in confronting Issue 2, we seek turbulent solutions 
that possess ft profiles with fast equators and strong 
latitudinal contrasts Aft, and emerging heat fluxes 
that vary little with latitude. To achieve this we 
have considered two paths in parameter space that 
yield more turbulent solutions by either varying the 
Prandtl number or keeping it fixed, while maintain- 
ing the same rotational constraint as measured by a 
convective Rossby number. 

We describe briefly in §2 the ASH code and the set 
of parameters used for the simulations studied here. 
In §3 we discuss the properties of rotating turbulent 
convection and the resulting differential rotation and 
the meridional circulation for the five cases A, AB, 
B, C and D. In §4 we analyze the transport of angu- 
lar momentum by several processes and the influence 
of baroclinic effects in establishing the mean flows. 
In §5 we reflect on the significance of our findings, 
and especially in terms of dealing with the two issues 
raised by the prior simulations with ASH. 

2. FORMULATING THE MODEL 

Our numerical models are intended to be a faithful 
if highly simplified descriptions of the solar convec- 
tion zone. In brief overview, solar values are taken 
for the heat flux, rotation rate, mass and radius, and 
a perfect gas is assumed since the upper boundary of 
the shell lies below the H and He ionization zones; 
contact is made with a real solar structure model 
for the radial stratification being considered. The 
computational domain extends from about 0.72i? Q 
to O.96i?0, where Rq is solar radius, with such shells 
having an overall density contrast in radius of about 
25, and as a concequence compressibility effects are 
substantial. Thus we are concerned only with the 
central portion of the convection zone, dealing neither 
with the penetrative convection below that zone, nor 
the two shear layers present at the top and bottom 
of it. Given the computational resources available, 
we prefer to concentrate our effort on processes that 
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establish the primary differential rotation in the bulk 
of the convection zone, and in future studies will seek 
to incorporate the other regions. We have as well 
softened the effects of the very steep entropy gradi- 
ent close to the surface that would otherwise favor 
the driving of very small granular and mesogranular 
scales of convection, with these requiring a spatial 
resolution at least ten times greater than presently 
available. 

The ASH code solves the 3-D anelastic equations of 
motion in a rotating spherical shell geometry using a 
pseudo-spectral semi-implicit approach (Clune et al. 
1999). As discussed in detail in Miesch et al. (2000), 
these equations are fully nonlinear in velocity vari- 
ables and linearized in thermodynamic variables with 
respect to a spherically symmetric mean state having 
a density p, pressure P, temperature T and specific 
entropy S and perturbations about this mean state 
of p, P, T, S. The conservation of mass, momentum 
and energy (or entropy) in a rotating reference frame 
are thus expressed as 



V • (pv) = 0, 



(1) 



dv 



P ( 0j + (v ■ V)v + 2H x v 



-VP 



+p g -V -©-[VP-jog], (2) 



as 



pT— = V-[K r pcpV(T + T) + KpTV(S + S)] (3) 



pTw ■ V(5 + S) + 2pv 



1/3(V • v) 



where c p is the specific heat at constant pressure, 
v = (vrjVQjVc/)) is the local velocity in spherical geom- 
etry in the rotating frame of constant angular velocity 
ft Q , g the gravitational acceleration, n r the radiative 
diffusivity, and V the viscous stress tensor, with com- 
ponents 



-2pu[e ij -l/3(y-v)S ij ], 



(4) 



where is the strain rate tensor. Here v and k are 



effective eddy diffusivities for vorticity and entropy. 
To close the set of equations, the linearized relations 
for the thermodynamic fluctuations are 
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assuming the ideal gas law 



P = KpT, 



(5) 



(6) 



where 7Z is the gas constant. The bracketed term in 
equation (2), VP — pg, vanishes initially because the 



mean state begins in hydrostatic balance from a one- 
dimensional radial solar model (Brun, Turck-Chieze 
& Zahn 1999), but as the convection becomes estab- 
lished this term becomes nonzero through effects of 
turbulent pressure. It is essential to take into account 
effects of compressibility upon the convection, since 
the solar convection zone spans many density scale 
heights. To accommodate this, we use the anelastic 
approximation (Gough 1969) to filter out the sound 
waves and therefore permit bigger time steps for the 
temporal evolution. The latter is allowed since the 
CFL (Courant, Friedrichs & Lewy) numerical stabil- 
ity condition now applies to the smaller convective 
velocities rather than the sound speed c s . 

Due to the small solar molecular viscosity, direct 
numerical simulations (DNS) of the full scale range 
of motions present in stellar convection zones are cur- 
rently not feasible. We seek to resolve the largest 
scales of convective motion that we believe are the 
main drivers of the solar differential rotation, doing 
so within a large-eddy simulation (LES) formulation 
where v and k are assumed to be an effective eddy 
viscosity and eddy diffusivity that represent unre- 
solved subgrid-scale (SGS) processes, chosen to suit- 
ably truncate the nonlinear energy cascade. For sim- 
plicity, both are here taken to be functions of radius 
alone, and are chosen to scale as the inverse of mean 
density. Other forms that may be determined from 
the properties of the large-scale flows according to 
one of many prescriptions (e.g. Lesieur 1997, Canuto 
1999) will be considered in the future. We have also 
introduced an unresolved enthalpy flux proportional 
to the mean entropy gradient in equation (3) in or- 
der to account for transport by small-scale convec- 
tive structures near the top of our domain (Miesch 
et al. 2000). We utilize the same radial profile for 
that mean eddy diffusivity in our five cases in order 
to minimize the impact of our SGS treatment on the 
main properties of our solutions. We emphasize that 
currently tractable simulations are still many decades 
away in parameter space from the intensely turbulent 
conditions encountered in the sun, and thus these 
large-eddy simulations must be viewed as training 
tools for developing our dynamical intuition of what 
might be proceeding within the solar convection zone. 

Within the ASH code, the mass flux is imposed to 
be divergence-free by using poloidal W and toroidal 
Z functions. The thermodynamic variables P and S, 
and W and Z, are expanded in spherical harmonics 
Yf l (6,(f)) to resolve their horizontal structures and 
in Chebyshev polynomials T n (r) to resolve their ra- 
dial structures. This approach has the distinct ad- 
vantage that the spatial resolution is uniform every- 
where on a sphere when a complete set of spherical 
harmonics is used in degree t (retaining all azimuthal 
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orders m). We expand up to degree £ = £ ma x (de- 
pending on the number of latitudinal mesh points 
Ng, e.g. l max = (2N e - l)/3), utilize as longitudi- 
nal mesh points = 2Ng, and employ N r colloca- 
tion points in projecting upon the Chebyshev poly- 
nomials. In this study the highest resolution used 
has £ m ax = 340 and N r = 193. The time evolution 
is carried out using an implicit, second-order Crank- 
Nicholson scheme for the linear terms and an explicit, 
second-order Adams-Bashforth scheme for the advec- 
tive and Coriolis terms. 

Within ASH, all spectral transformations are ap- 
plied to data local to each processor, with inter- 
processor transposes performed when necessary to ar- 
range for the transformation dimension to be local. 
The triangular truncation in spectral space precludes 
any simple distribution of the data and workload 
among the nodes. For very large problems, the Leg- 
endre transformations dominate the workload and, 
as a result, great care has been taken to optimize 
their performance on cache-based architectures. Ar- 
rays and loops have been structured to operate on 
blocks which minimize cache misses. The ASH code 
is extremely flexible and has demonstrated excellent 
scalability on massively parallel supercomputers such 
as the Cray T3E, IBM SP-3 and Origin 2000. 

As boundary conditions, we impose impenetrable 
and stress-free conditions for the velocity field and 
constant flux (i.e constant entropy gradient) both at 
the inner and outer boundaries. We seek solutions 
with an emerging flux at the top that is invariant 
with latitude (Issue 2). As initial conditions, we have 
started some simulations (cases A, B) from quiescent 
conditions of uniform rotation, and others (cases AB, 
C, D) from evolved solutions in which we modify cer- 
tain diffusivities. This leads to changes in the effec- 
tive Rayleigh number R a , the Prandtl number P r , the 
Peclet number P e , the Reynolds number R e and the 
Taylor number T a , while keeping constant the con- 
vective Rossby number R oc , all of which are defined 
in Table 1. We also summarize there the parameters 
of the five simulation cases. 

3. PROPERTIES OF TURBULENT COMPRESSIBLE 
CONVECTION 

We have conducted five simulations involving in- 
creasingly nonlinear flows that are achieved by reduc- 
ing the viscous and entropy diffusivities in the manner 
outlined in Table 1. We have followed two paths in 
parameter space in obtaining more complex convec- 
tive flows. On Path 1 in going from case A to C via 
B, we incrementally decreased the eddy viscosity v 
while keeping the eddy diffusivity n constant, thereby 
reducing the Prandtl number P r by a factor 8. In par- 
ticular, the laminar case A has P r of unity; reducing 
the viscosity by a factor of 4 leads to the mildly tur- 



bulent case B with P r = 0.25, or by a factor of 8 leads 
to the more turbulent case C with P r = 0.125. This 
serves to increase the Reynolds number R e while only 
mildly increasing the Peclet number P e . Path 2 kept 
the Prandtl number fixed at P r = 0.25 as the com- 
plexity of the flows was increased by reducing both 
diffusivities. Starting from case AB, we go to case B 
by decreasing both diffusivities v and n by a factor of 
2, and then to our most turbulent case D by further 
reducing both by a factor of 2 relative to case B. This 
Path 2 in going from case AB to D via B results in 
both R e and P e increasing comparably. All our mod- 
els possess a convective Rossby number R oc of order 
2/3, thus maintaining a strong rotational constraint 
on the convection. 

As we shall describe in some detail, the resulting 
vigorous convection influenced by rotation in all these 
cases is intricate and richly time dependent, much as 
found in Miesch et al. (2000) and Elliott et al. (2000). 
It is characterized by networks of strong downflow 
at the periphery of the convection cells, and weaker 
upflows in their middle, both of which are a conse- 
quence of the effects of compressibility as we consider 
flows that can span multiple density scale heights in 
the vertical. Indeed, we consistently observe that the 
downflows are able to extend over the full depth of the 
unstable layer, appearing as twisted sheets of down- 
flow near the top and more distinctive plumes deeper 
in the layer. These downflow networks essentially 
represent coherent structures amidst the turbulence, 
and they are found to have a most significant role 
in the nonlinear transport of angular momentum by 
yielding correlations between different velocity com- 
ponents that form Reynolds stress terms. We find 
that the convection in all cases studied here is able 
to redistribute angular momentum in such a manner 
that substantial differential rotation profiles are es- 
tablished, the properties of which are the major focus 
of this work. 

3.1. Complex evolution of convective patterns 

The time dependence in our most turbulent simula- 
tion (case D) is shown in Figure [2] which displays two 
sequence of images of the radial velocity on spheri- 
cal surfaces over the course of one full rotation. The 
upper sequence with views near the top of the layer 
involves simpler downflow networks (shown in darker 
tones) that are easier to intercompare from frame to 
frame, whereas the lower ones with views in the mid- 
dle of the convecting layer are more difficult to track 
because of increased complexity of the patterns in the 
more turbulent flows there. The vantage point is in 
the uniformly rotating frame used in our modelling, 
and some of the pattern evolution results from the 
prograde zonal flows at low latitudes and retrograde 
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Fig. 2. — Evolution of the convection throughout one solar rotation, showing the radial velocity of case D near the top and at 
the middle of the domain. The time interval between each successive image is about 7 days. Features 1 and 2 exhibit the merging 
of the downflow lanes, feature 3 the shearing action of the differential rotation present in the shell and feature 4 the appearance 
and deformation of a convective cell at higher latitudes. 



ones at high latitudes associated with the differen- 
tial rotation relative to this frame. There is further 
melding and shearing of particular downflow lanes as 
the convection cells evolve over a broad range of time 
scales, some of which are comparable to the rota- 
tion period. This is particularly evident in some of 
the downflow structures identified near the equato- 
rial region in the upper sequence, with features la- 
beled 1 and 2 illustrating the merging of two down- 
flow lanes, and feature 3 the typical distortion of a 
lane which also involves both a site of cyclonic swirl 
in the northern hemisphere and another that is ap- 
propriately anticyclonic in the southern hemisphere. 
The behavior at higher latitudes that involves retro- 
grade displacement of the downflow networks is some- 
what more intricate, partly because the convection 
cells are of smaller scale and exhibit the frequent for- 
mation of new downflow lanes (as in feature 4) that 
can serve to cleave existing cells. Figure § empha- 
sizes that the overall pattern of these global cells is 
sufficiently modified during the course of one rota- 
tion period so that it would be difficult to identify 
particular structures (relative to our uniformly ro- 
tating vantage point) when viewed in a subsequent 
rotation. This would suggest that giant cells pos- 
sibly present within the solar convection zone may 
also loose their identity from one Carrington rotation 
to the next. This comes about both because of ad- 
vection and distortion of the cells by the mean zonal 
flows associated with the differential rotation (here at 
the equator leading to relative angular displacements 



in longitude of about 70° over one rotation period), 
and because of fairly rapid evolution and some prop- 
agation in their individual downflow patterns. 

3.2. Downflow networks and variation with depth 

The convective structures as delineated by the 
downflow networks show distinctive changes as the 
level of turbulence is increased in going from case A to 
case D. Figure |3] provides an overview of radial veloc- 
ity snapshots in our five simulations at three depths 
(near the top, middle and bottom), accompanied by 
the fluctuating temperature fields at mid-depth. The 
upper surface in all our cases involves a connected 
network of downflows surrounding broad upflows, but 
such smoothness can disguise far more turbulent flows 
below. The seemingly cellular motions near the sur- 
face result from the expansion of fluid elements rising 
through the rapidly decreasing density stratification 
near the upper boundary, aided also by our increas- 
ing viscous and thermal diffusivities there. As viewed 
near the top, the tendency of the convection in our 
laminar case A to be organized into 'banana cells' 
nearly aligned with the rotation axis at low latitudes 
is progressively disrupted by increasing the level of 
complexity in going in turn to cases AB, B, C and 
D. There is still some semblance of north-south align- 
ment in the downflows even in our most turbulent 
case D, but the latitudinal span of this alignment 
is confined to a narrow interval around the equator. 
Clearly the downflow lanes become more wiggly and 
exhibit more pronounced vortical features and curva- 
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Fig. 3. — Convective patterns for the five cases A, AB, B, C and D as increasingly turbulent flows are attained. The radial 
velocity shapshots are shown at three different depths (0.95, 0.84, 0.73 Rq). Downflows are represented in purple-dark tones and 
upflows in orange-bright tones, with dynamic ranges indicated. The dotted circle is positioned at radius R.0, and the equator is 
indicated by the dashed curve. The convective structures become more complex in this progression of cases, with the banana-like 
convective cells giving way to stronger and more frequent vortex sites. The strongest downflow lanes extend over the full depth 
range. The fluctuating temperature fields at mid depth are shown on the right, emphasizing that downflows are relatively cool 
and that the polar regions are on average warm. 
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ture in this sequence of cases. As well, the downflow 
networks involve more frequent branching points and 
smaller horizontal scales for the convective patterns, 
especially at higher latitudes. Given the three simul- 
taneous views of the radial velocity, one can clearly 
identify downflow lanes near the top in all our cases 
that turn into distinctive plumes at greater depths, 
showing that organized flows extend over multiple 
scale heights. Indeed, the strongest downflows oc- 
cur at the interstices of the upper network and are 
able to pierce through the interior turbulence, thus 
spanning the full depth range of the domain. 

The plumes in the more turbulent cases C and 
D represent coherent structures that are embedded 
within less ordered flows that surround them. They 
are able to maintain their identity, though with some 
distortion and mobility, over significant intervals of 
time. Although these downflowing plumes are pri- 
marily directed radially inward, they show some tilt 
both toward the rotation axis and out of the merid- 
ional plane. This yields correlated velocity compo- 
nents and thus Reynolds stresses that are a key in- 
gredient in the redistribution of angular momentum 
within the shell. Such tilting away from the local ra- 
dial direction in coherent downflows has been seen 
in high-resolution local /-plane simulations of rotat- 
ing compressible convection (Brummell et al. 1998), 
and their presence has a dominant role in establishing 
the mean zonal and meridional flows. We also refer 
to Rieutord & Zahn (1995) and Zahn (2000) for an 
analytical study of the transport properties and cor- 
relations present in such strong vortex structures and 
on their potential dynamical role in the solar convec- 
tion zone. 

The strong downflows shown in Figure || accentu- 
ate the asymmetries that are characteristic of com- 
pressible convection, with typical peak amplitudes 
in these downflows at mid-layer being as much as 
twofold greater than that in the upflows. As might 
be expected, the overall rms radial velocities listed 
in Table 2 increase with complexity in the flow fields 
in going from case A to D. The asymmetries be- 
tween upflows and downflows have the consequence 
that the kinetic energy flux in such compressible con- 
vection is directed radially inward, in contrast to the 
enthalpy and radiative fluxes that are directed out- 
ward in transporting the solar flux (see Figure 10 a). 

That enthalpy flux involves correlations between 
radial velocities and temperature fluctuations, and 
these are evidently strong as seen when inspecting the 
temperature and velocity fields shown at mid layer 
in Figure [3|. Buoyancy driving within our thermal 
convection involves downflows that are cooler and 
thus denser and upflows that are warmer and lighter 
than the mean; there are systematic asymmetries in 



those temperature fluctuations, much as in the ra- 
dial velocities. Further, in comparing the tempera- 
ture maps with those of radial velocity in the mid- 
dle of the layer, some of the temperature patterns 
are evidently smoother, which is a consequence of 
the greater thermal diffusivities than viscosities for 
cases with Prandtl numbers less than unity. A strik- 
ing property shared by all these temperature fields is 
that the polar regions are consistently warmer than 
the lower latitudes, a feature that we will find to be 
consistent with a fast or prograde equatorial rotation. 

3.3. Driving strong differential rotation 

The differential rotation profiles with radius and 
latitude that result from the angular momentum re- 
distribution by the vigorous convection in our five 
simulations are presented in Figure In order to sim- 
plify comparison of our results with deductions drawn 
from helioseismology (Fig.Q) , we have converted our 
mean longitudinal velocities {u (with " denoting aver- 
aging in longitude and time) into a sidereal angular 
velocity f2 with radius and latitude, and note that 
our reference frame rotation rate Q /2ir is 414 nHz 
(or a period of 28 days). The angular velocity in all 
our cases exhibits substantial variations in time, and 
thus long time averages must be formed to deduce 
the time mean profiles of Q shown in Figure The 
layout of the five cases in Figure |I| reflects the two 
paths we have taken in increasing the complexity or 
turbulence level in the solutions: Path 1 in going from 
A to C via B while decreasing the Prandtl number 
takes us from upper left to lower right, and Path 2 
in going from case AB to D via B while keeping the 
Prandtl number fixed at P r =l/4 takes us from upper 
right to lower left. Complexity in the convective flows 
increases in going down the page. 

All five simulations yield angular velocity f2 profiles 
that involve fast (prograde) equatorial regions and 
slow (retrograde) high latitude regions. The variation 
of with radius and latitude may be best judged in 
the color contour plots in Figure f| which are scaled in- 
dependently for each of the cases; the reference frame 
rate is also indicated. The immediate polar regions 
are omitted in these plots because it is difficult to 
obtain stable mean f2 values at very high latitudes 
since the averaging domain there becomes very small 
whereas the temporal fluctuations in the flows remain 
substantial. The contour plots reveal that there are 
some differences in the f2 realized in the northern and 
southern hemispheres, though such symmetry break- 
ing is modest and probably will diminish with longer 
averaging. The convection itself is not symmetric 
about the equator, and thus the mean zonal flows 
that accompany such convection, and which are man- 
ifest as differential rotation, can be expected to have 
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Fig. 4. — Temporal and longitudinal averages of the angular velocity profiles achieved in cases A, AB, B, C and D, formed 
over intervals in turn of 295, 275, 275, 175 and 35 days. The contour plots for Q/2n on the left of each panel are independently 
scaled, whereas the radial cuts at the indicated latitudes share the same scaling to accentuate the overall behavior of the five 
cases. The crossed layout of the five cases emphasizes the two different paths followed to reach more turbulent states, mainly by 
lowering P r on Path 1 (A — * B — » C), and by lowering diffusivities while keeping P r constant on Path 2 (AB — > B — > D). All 
cases exhibit a prograde equatorial rotation and a strong contrast AQ from equator to pole. Case AB possesses a high latitude 
region of particularly slow rotation. 
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variations in the two hemispheres. In cases B, C and 
D, there is some alignment of the f2 contours at the 
lower latitudes with the rotation axis, thus showing 
a tendency for fi to be somewhat constant on cylin- 
ders. Further, in these cases almost all the decrease 
in f2 with latitude occurs in going from the equator to 
about 45°, or thus is confined to the region outside 
the tangent cylinder to the inner boundary (which 
intersects the outer boundary in our shell configu- 
ration at 42°). In contrast, cases A and AB show 
far less alignment of O contours with cylinders at the 
lower latitudes, and at mid-latitudes the contours are 
nearly aligned with radial lines, more in the spirit of 
the hehoseismic inferences. 

Case AB in Figure |4| is unique in having the mono- 
tonic decrease of with latitude continue onward to 
high latitudes, which is also the trend deduced from 
hehoseismic measurements. Thus Issue 1 concerned 
with achieving a consistently decreasing $7 at high 
latitudes is resolved with case AB. This is significant 
in showing that such behavior can be realized in our 
modelling of convection in deep shells, though it is not 
a common property in our other cases. It would be 
most desirable to understand how such high-latitude 
variation in f2 is achieved in case AB, and we will 
address this in §4. 

The accompanying radial cuts of at six fixed 
latitudes in Figure |] permit us to readily quantify 
the £1 contrasts with latitude achieved in these so- 
lutions, and to judge the functional variation with 
radius in each case. We use a common scaling for all 
these line plots to make intercomparison between the 
cases most convenient; the radial cuts for f2 have been 
averaged between the north and south hemispheres. 
Near the top of the convection zone at radius 0.96i?, 
the laminar convection in case A produces a differ- 
ential rotation with a contrast in angular velocity, or 
A$7/27r, of about 50 nHz between the equator and 
60°, or 12% relative to the frame rotation rate (also 
quoted in Table 2). Continuing on Path 1 in param- 
eter space to the more turbulent cases B and C, we 
find that the latitudinal contrast in angular velocity 
has increased substantially, becoming 115 nHz and 
125 nHz in the two cases. These correspond in turn to 
a 28% and a 30% variation of the rotation rate. These 
values are of interest since the hehoseismic inferences 
(Thompson et al. 1996, Schou et al. 1998, Howe et al. 
2000b) have a contrast of about 92 nHz at a similar 
depth between the equator and 60°, or a variation of 
about 22% in rotation rate, which further increases 
to about 32% in going to 75°. The pronounced differ- 
ential rotation in cases B and C is accompanied by 
the O profiles becoming somewhat more aligned with 
the rotation axis, resulting in steeper slopes in the 
radial cuts at low and mid latitudes. These two tur- 



bulent cases achieve their larger A£l by both faster 
equatorial rotation rates and slower rates at higher 
latitudes. Thus Path 1 has been able to resolve Issue 
2, concerned with retaining a strong contrast Af2 and 
a fast equator, as the solution becomes more complex 
and turbulent. 

Turning to Path 2 in parameter space, case AB 
shows a contrast of about 135 nHz between the equa- 
tor and 60°, or a 33% variation of rotation rate, which 
further increases to about 160 nHz or 39% in going 
to 75°. The pivotal case B has a somewhat reduced 
contrast AO/27T of 115 nHz or 28% variation between 
equator and 60° , with little further variation at higher 
latitudes. The most turbulent case D has a AO/2-7T of 
about 105 nHz or a 25% variation between the equa- 
tor and 60° . Thus Path 2 leads to a slight reduction in 
AQ with increasing complexity, unlike the behavior 
of Path 1. However, even this path yields a turbu- 
lent solution case D whose Ail is still close to the 
hehoseismic contrast, thus largely resolving Issue 2. 
This is reemphasized in Figure |H| which summarizes 
the variation of Ail with P r for our five cases. 
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Fig. 5. — Parameter space diagram for relative latitudinal 
angular velocity contrast Afi/fi as a function of the Prandtl 
number P r for the five cases. The two paths toward higher lev- 
els of turbulence either reduce P r (A — > B — > C), or maintain 
a constant P r {AB — > B — > D). 

Most of our cases possess overall latitudinal con- 
trasts in f2 that are in the realm of solar values de- 
duced from inversion of hehoseismic data, yet case 
AB stands out in having the systematic decrease of 
f2 with latitude extending almost to the poles, which 
appears to be another distinguishing feature of the 
actual solar $7 profiles. Further, case AB displays 
little radial variation in Q at intermediate and high 
latitudes (from say 45° onward) as the angular ve- 
locity continues to decrease poleward. Such behavior 
is most interesting, and it is necessary to understand 
just which convective properties within case AB allow 
it to come into reasonable contact with the hehoseis- 
mic profiles for £1 deduced in the bulk of the solar 
convection zone. 
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The Q profiles in Figure || have been formed from 
temporal averages spanning multiple rotation peri- 
ods as indicated. It is appropriate to consider if 
these represent truly 'spun-up' solutions in a statis- 
tical sense, and further whether several distinctive O 
profiles could be achieved for the same control pa- 
rameters. Both issues may be intertwined, for the 
rate of approach to equilibration can be influenced 
by the attraction characteristics of that differential 
rotation state, and of course by the amplitude of the 
fluxes available to redistribute angular momentum to 
achieve that state (see §4.1). This overall dynamical 
system of turbulent convection is sufficiently complex 
that we are uncertain whether there may exist mul- 
tiple basins of attraction leading to different classes 
of differential rotation. For instance, is the behav- 
ior of case AB with noticeably slow rotation at high 
latitudes an example of one class of behavior, and 
our other cases that of another family? Could such 
families overlap in parameter space, or are there just 
gradual variations in behavior in Q with changes in 
the parameters? We have so far sought to address 
some of these questions by perturbing the evolving 
solutions to see if they might flip to another state, 
but they have not done so. We plan to examine such 
issues of solution uniqueness in our follow on studies 
in which we seek to extend the slow-pole characteris- 
tics of case AB to other parameter settings involving 
more complex convection. 

As to the relative maturity of the spun-up states 
shown in Figure ||, these vary from case to case due 
to the rapidly increasing computational expense in 
dealing with the finer spatial resolution required by 
the more complex simulations. Cases AB and C were 
both started from case B that had already been run 
for over 4000 days of elapsed simulation time (or a 
nominal 143 rotation periods involving about 28 days 
each). At this point case B appeared to be statisti- 
cally stationary in terms of the kinetic energy asso- 
ciated with the differential rotation, though it like 
most of the other simulations exhibit small fluctu- 
ations in Vt profiles determined from single-rotation 
averages, especially at the higher latitudes. Case AB 
was evolved for about 2300 days (82 rotation peri- 
ods), and we illustrate in Figure ga a succession of 
Q profiles with latitude sampling the last 600 days in 
the simulation. The solid curve there is an average 
formed over 10 rotations (consistent with the con- 
tour plot in Figure , and we see that the individual 
rotation averages being sampled form a narrow en- 
velope around it. There is evidently some symmetry 
breaking between the two hemispheres. We believe 
that the differential rotation for case AB is now an 
effectively stationary state (as is also confirmed in 
studying the angular momentum flux balance in Fig- 



ure [ll]). The more turbulent case C was evolved for 
about 500 days (18 rotations) after being initiated 
from case B, and a set of its angular velocity profiles 
are shown sampling the last 300 days in Figure |6^. 
We are less certain of its stationarity, but we could 
not detect any systematic trends in the evolution of 
its differential rotation over the last 10 rotations. We 
saw no evidence of a slow pole developing, but that 
may well require more extended computations than 
could be presently arranged. Figure |6| serves to em- 
phasize that the angular velocity even in the sun may 
be expected to vary somewhat from one rotation to 
another, with the samplings here providing a sense of 
the amplitude of those changes. 
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Fig. 6. — Succession of time-averaged Q profiles with lati- 
tude at r = O.96J?0. (a) For case AB, a numbered sequence 
of single-rotation averages spanning an interval of 600 days in 
the late evolution of the system, with the bold curve 4 denot- 
ing an average over the last 275 days in the simulation, (b) 
For case C, dealing with samples in a 300 day interval, and the 
bold curve 4 representing an average over the final 175 days. 
The variations are representative of small changes in the dif- 
ferential rotation that accompany changes in the convection 
once a mature statistical state has been achieved. 

3.4. Meridional circulation patterns 

The time-averaged meridional circulations that ac- 
company the vigorous convection in the five cases 
are shown in Figure [7[ The typical amplitudes in 
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Fig. 7. — Temporal and longitudinal averages of the meridional flows achieved in the cases A, AB, B, C and D, deduced from 
sampling in turn 295, 275, 275, 175 and 35 days. Shown are random streaklines whose length is proportional to flow speed, with 
arrowheads indicating flow sense. The typical speeds in these meridional circulations are about 20 ms _1 . For all the cases strong 
poleward cells are present near the surface at low latitudes as well as return flows at mid depth. 



these large-scale circulations are about 20 ms , and 
thus comparable to the values deduced from local 
domain helioseismic probing of the uppermost con- 
vection zone based either on time-distance methods 
(e.g. Giles, Duvall & Scherrer 1998) or ring-diagram 
analyses (Schou & Bogart 1998, Haber et al. 1998). 
There is little change in meridional circulation ampli- 
tudes as we increase the level of turbulence in going 
from case A to D. However, multi-cell structures in 
these circulations become more intricate with the in- 
creased complexity of the convection. At lower lati- 
tudes the circulation in both hemispheres is poleward 
near the top of the domain, with return flows at var- 
ious depths. All cases display multiple cells with ra- 
dius and latitude, and never only one big meridional 
cell as often used in mean field models dealing with 
differential rotation (e.g. Rekowski & Riidiger 1998, 
Durney 2000) or with Babcock-Leighton dynamos 
(e.g. Choudhuri, Schiissler & Dikpati 1995, Dikpati 
& Charbonneau 1999). The resulting axisymmet- 
ric meridional circulation is maintained by Coriolis 
forces acting on the mean zonal flows that appear 
as the differential rotation, by buoyancy forces, by 
Reynolds stresses, and by pressure gradients. Given 
these competing processes, it is not self evident as 
to what pattern of circulation cells should result, nor 
how many should be present in depth or latitude. 
Our five simulations have shown that there is some 
variety in the meridional circulations achieved, all of 
which involve multi-celled structures. Since the ki- 
netic energy in the meridional circulation (MCKE) is 
typically about two orders of magnitude smaller than 
in the differential rotation (DRKE), as we will detail 
in §3.5 and in Table 2, small variations in the differ- 



ential rotation can yield substantial changes in the 
circulations. This is likewise true of the time-varying 
Reynolds stresses from the evolving convection which 
again has a kinetic energy (CKE) much larger than 
that of the meridional circulations. This may explain 
the complex time dependence realized by the merid- 
ional flows, and the need to use long time averages in 
defining their mean properties. 

Another rendition of the time-averaged meridional 
circulations achieved in cases AB and C is shown 
in Figure || using a streamfunction $ based on the 
zonally-averaged mass flux [as in Miesch et al. 2000, 
equation (7)]. In case AB (Fig. ^|a) there are two cir- 
culation cells positioned above each other in radius 
at low latitudes. The stronger upper one (solid con- 
tours representing counterclockwise circulation) in- 
volves poleward flow that extends from the equator 
to about 30° latitude near the top of the domain in 
the northern hemisphere. The southern hemisphere 
has likewise poleward flow near the top at low lati- 
tudes, with ascending motions again present from the 
equator to about 20° latitude. At latitudes greater 
than about 30° the relatively weak flow near the top 
is mainly equatorward in both hemispheres, but ex- 
hibits fluctuations. 

A quantitative measure of this for case AB is pro- 
vided in Figure that displays the mean velocity 
component v$ with latitude at two depths near the 
top of the domain. The poleward flow in both hemi- 
spheres peaks at about 20° latitude and then de- 
creases rapidly, changing to weak equatorward flow 
above 30° which attains about one-third that peak 
amplitude. 
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Fig. 8. — Streamlines of the mean axisymmetric merid- 
ional circulation achieved in (a) case AB averaged over 275 
days, and in (b) case C averaged over 175 days. Solid contours 
denote counterclockwise circulation (and dashed clockwise), 
equally spaced in value. In case AB, two circulation cells are 
present with radius at low latitudes, and only weak circula- 
tions at latitudes above 30°. Case C possesses three cells at 
low latitudes, with the deepest extending prominently to high 
latitudes. 

Turning to case C in Figure ||&, it exhibits three 
circulation cells positioned radially at low latitudes, 
with the outermost again yielding poleward flow at 
the top of the domain that extends to about 35° in 
latitude. At higher latitudes the mean meridional 
flow is again equatorward near the top, attaining a 
peak amplitude for vg (detailed in Figure that is 
comparable to the poleward one from the low lati- 
tudes, unlike in case AB. Of the three meridional 
cells at low latitudes in Case C, much as for model 
TUR in Miesch et al. 2000 (cf. Fig. 16a), the deep- 
est cell involves a strong counterclockwise circulation 
that extends to high latitudes, yielding a submerged 
poleward flow there that lies below the equatorward 
flow at the top of the domain. Such behavior involv- 
ing a third deep circulation cell that extends to high 
latitudes is also seen in cases B and D. Such a strong 
third cell appears to be of significance in the contin- 
uing net poleward transport of angular momentum 
by the meridional circulations (cf. §4.1, Fig. in 
all these cases at latitudes above about 30°. This is 
not realized in case AB, and may contribute to its 
slow pole behavior. It is encouraging that we have 
poleward circulations in the upper regions of the sim- 
ulations, which is in accord with the general sense of 
the mean flows near the surface being deduced from 
local helioseismology, though two cell behavior with 
latitude has been detected recently only in the north- 
ern hemisphere near the peak of solar activity (Haber 
et al. 2000). Such symmetry breaking in the two so- 
lar hemispheres is an interesting property, and one 
that is also occasionally realized in our simulations 
as the convection patterns evolve. The helioseismic 
probing with ring diagram methods and explicit in- 



versions is able to sense the meridional circulations 
only fairly close to the solar surface, typically extend- 
ing to depths of about 20 Mm or to radius 0.97 Rq, 
whereas our simulations have their upper boundary 
slightly below this level at 0.96 Rq. 
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Fig. 9. — Mean velocity component ve with latitude at the 
two depths r = 0.96i?o (solid) and O.94i?0 (dash dot), showing 
in (a) case AB and in (b) case C. Positive values correspond 
to flow directed from north (positive latitudes) to south (neg- 
ative latitudes). At low latitudes the flows are poleward in 
both hemispheres, but whereas case C exhibits fairly strong 
equatorward flow at latitudes above 35°, case AB possesses 
much weaker circulations there. 

Thus we must be cautious in interpreting similar 
behavior in the meridional circulations since our mod- 
els and the ring diagram analysis do not explicitly 
overlap in radius. Helioseismic assessments based 
on time-distance methods (Giles 1999, Chou & Dai 
2001) and annular rings centered on the poles (Braun 
& Fan 1998) report detecting effects attributable to 
meridional circulations with a mainly poleward sense 
to depths corresponding to 0.90 Rq or even 0.85 Rq. 
Such results are most interesting, but considerable 
further work on inversions would be required to pro- 
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vide detailed profiles of the circulations with depth. 
As these mappings become available, they may be 
able to confirm or refute the multi-cell radial struc- 
ture of meridional circulation (Fig. |7|) typically real- 
ized in our simulations. 

3.5. Energetics of the convection and the mean flows 

The overall energetics within these shells of rotat- 
ing convection have some interesting properties, in 
addition to the mean zonal and meridional flows that 
coexist with the complex convective motions. The 
convection is responsible for transporting outward the 
solar flux emerging from the deep interior. We should 
recall, as discussed in detail in Miesch et al. (2000), 
that the radial flux balance in these convective shells 
involves four dominant contributors, namely the en- 
thalpy or convective flux F e , the radiative flux F r , 
the kinetic energy flux F^, and finally the unresolved 
eddy flux F u , which add up to form the total flux 
Ft- Figure ICa shows the flux balance with radius 



achieved in our most turbulent case D as averaged 
over horizontal surfaces and converted to luminosi- 
ties. The radiative flux becomes significant deep in 
the layer due to the steady increase of radiative con- 
ductivity with depth, and indeed by construction it 
suffices to carry all the imposed flux through the 
lower boundary of our domain where the radial ve- 
locities and thus the convective flux vanishes. A sim- 
ilar role near the top of the layer is played by the 
sub-grid-scale turbulence that yields F u , which being 
proportional to a specified eddy diffusivity function 
k and the mean radial gradient of entropy, suffices to 
carry the total flux through the upper boundary and 
prevents the entropy gradient there from becoming 
too superadiabatic compared to the scales of convec- 
tion that we are prepared to resolve spatially. Over 
most of the interior of the shell, the strong correla- 
tions between radial velocities and temperature fluc- 
tuations yield the enthalpy flux F e that transports 
upward almost all of the imposed flux, and this peaks 
near the middle of the layer. The kinetic energy 
flux Fk works against the others by being directed 
downward, a result of the fast downflow sheets and 
plumes achieved by effects of compressibility (Hurl- 
burt, Toomre & Massaguer 1986). These general 
properties are shared by our five cases, all of which 
have achieved good overall flux balance with radius, 
as can be assessed by examining Ft. 

Figure [lO| b presents the kinetic energy spectra with 
azimuthal wavenumber m at three depths, and aver- 
aged in time, as realized in the case D simulation. 
The spectra are fairly broad, with a plateau of power 
extending up to about m=30 corresponding to some 
of the most vigorously driven scales, and then a rapid 
decrease involving about 5 orders of magnitude to the 



highest wavenumber of 340. The decrease is more 
rapid for the spectra formed near the top of the shell. 
These spectra suggest that the flows are well resolved, 
with a reasonable scale separation between the dom- 
inant energy input range and the wide interval over 
which dissipation functions. We cannot readily iden- 
tify a clear inertial subrange, though for reference 
we include some power laws. We also refer to Hath- 
away et al. (1996, 2000) for a discussion of recent 
observational inferences about the solar kinetic en- 
ergy spectrum which does not seem to indicate any 
clear scaling law. 




0.0 

-0.2 L_ 



0.75 0.80 0.85 0.90 0.95 

r/R 




10 100 
Azimuthal wavenumber m 

Fig. 10. — (a) Radial transport of energy in case D achieved 
by the fluxes F r , F e , Fk and F u , and their total Ft, all nor- 
malized by the solar luminosity, (b) Time averaged azimuthal 
wavenumber spectra of the kinetic energy on different radial 
surfaces (solid curve r/i?0=O.95; dashed, 0.84; dotted, 0.73) 
for case D with l ma ,x = 340. The results have been averaged 
over 35 days. Superimposed are the power laws k~°^ , k~ 2 
and k~ 3 ; no clear inertial subrange is identifiable. 

Table 2 summarizes various rms velocities that 
characterize our five simulations as sampled in the 
middle of the layer where the enthapy flux also peaks. 
The rms radial velocity v r increases monotonically in 
going through the sequence of cases A, AB, B, C 
to D. The associated rms Reynolds number R e in 
Table 1 increases also (though part of this is due to 
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changes in the diffusivities), varying by a factor of 
about 15 from our laminar to most turbulent solu- 
tions. The rms longitudinal velocity has the great- 
est amplitude in all the cases. However a removal of 
the mean zonal flow component responsible for the 
differential rotation yields v'^. Comparison of this 
with the radial and latitudinal rms velocities reveals 
that all possess very comparable amplitudes, suggest- 
ing fairly isotropic convective motions near the mid- 
plane. Table 2 also assesses the volume and time av- 
eraged total kinetic energy (KE), that associated with 
the differential rotation (DRKE), with the meridional 
circulation (MCKE), and with the convection itself 
(CKE). In all of our solutions the DRKE and CKE 
are comparable, and the MCKE is much smaller. Ta- 
ble 2 reveals a most interesting contrast in behavior 
for the two paths. Following Path 1 (involving A, B 
and C, with decreasing P r ), we find that KE in the 
solutions increases steadily with increasing flow com- 
plexity. This would be expected since the buoyancy 
driving has strengthened relative to the dissipative 
mechanisms as measured by the increasing Rayleigh 
number R a (Table 1). Path 2 (involving AB, B, and 
D, with P r kept fixed at 0.25) is quite different as 
R a increases: here the total kinetic energy KE re- 
mains nearly constant. A consequence is that with 
increasing complexity and increasing CKE along Path 
2, the DRKE must in turn decrease, and AO becomes 
smaller. This striking property of achieving a nearly 
constant KE along Path 2 (where both R e and P e in- 
crease comparably) is a remarkable feature of this in- 
tricate rotating system that is currently unexplained. 

Our solutions typically exhibit small differences in 
behavior in the two hemispheres, as can be detected 
in the time-averaged Q contours shown in Figure ||, 
and in the associated latitudinal cuts at fixed ra- 
dius displayed in Figure |6] for cases AB and C . The 
meridional circulations likewise show some symme- 
try breaking in their response between the northern 
and southern hemispheres in Figure 0, which is fur- 
ther quantified for cases AB and C in showing the 
meridional streamlines in Figure | and m examining 
the latitudinal variation of the mean velocity compo- 
nent vq in Figure |9|. A sense of these asymmetries 
can also be assessed by examining differences in the 
kinetic energy of differential rotation in the two hemi- 
spheres. For case AB, DRKE in the northern hemi- 
sphere is2.12xl0 6 erg cm 3 and 2.09 x 10 6 erg cm 3 
in the southern hemisphere, or a 1.6% difference. For 
case C, the corresponding values are 1.82 x 10 6 and 
1.76 x 10 6 , or 3.6%. We expect that such symmetry 
breaking is likely to evolve slowly, with neither hemi- 
sphere favored. We plan to study aspects of sym- 
metry breaking further with more extended simula- 
tions in the near future. Such efforts are inspired by 



the evolving meridional circulations and mean zonal 
flows being detected by helioseismology (Haber et al. 
2000, 2001), and the differing solar rotation rates in 
the two hemispheres deduced from tracking sunspots 
(Howard, Gilman & Gilman 1984). 

4. INTERPRETING THE DYNAMICS 

Our shells of rotating compressible convection are 
very complicated dynamical systems in terms of the 
nonlinear feedbacks and couplings that operate. It 
is difficult from first principles to predict or explain 
their overall behavior in terms of the differential rota- 
tion and meridional circulations that can be achieved 
and sustained as we sample different sites in param- 
eter space. The five simulations represent numerical 
experiments that seek to probe some of the families 
of responses within a highly simplified version of the 
solar convection zone. Although most of our approxi- 
mations here seem reasonable and necessary to yield a 
problem tractable to computational experiments, we 
do not fully know their impact and thus must draw 
our interpretations about the operation of the overall 
dynamics with considerable caution. The numerical 
solutions have the enormous advantage that we can 
interrogate them in detail to study various balances 
and fluxes, and these help to provide insights about 
the dynamical system. 

4.1. Redistributing the angular momentum 

Our choice of stress-free boundaries at the top and 
bottom of the computational domain has the advan- 
tage that no net torque is applied to our convective 
shells resulting in the conservation of the angular mo- 
mentum. We seek here to identify the main physical 
processes responsible for redistributing the angular 
momentum within our rotating convective shells, thus 
yielding the differential rotation seen in our five cases. 
We may assess the transport of angular momentum 
within these systems by considering the mean radial 
(J- r ) an d latitudinal (To) angular momentum fluxes. 
As discussed in Elliott et al. (2000), the (^-component 
of the momentum equation expressed in conservative 
form and averaged in time and longitude yields 

1 d(r^ r ) | 1 djsmdTe) _ (?) 
r 2 dr rsin^ 89 

involving the mean radial angular momentum flux 

T r = pr sin 9[-ur-^- (^-^j + v' r v'^ + v r (v ( p + Qor smO)} 

(8) 

and the mean latitudinal angular momentum flux 

_ „ . . r sin0 d ( Va, \ ~t~~i „ 

Tg = prsm0[-u — )+v g v (j) +v e (v ( f ) +n rsm 

r ov \ sin v J v 

(9) 
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In the above expressions for both fluxes, the first 
terms in each bracket are related to the angular mo- 
mentum flux due to viscous transport (which we de- 
note as T r y and Fey), the second term to the trans- 
port due to Reynolds stresses (J>,il and J 7 6,r) and 
the third term to the transport by the meridional 
circulation {T r ,M and Tq^m)- The Reynolds stresses 
above are associated with correlations of the velocity 

components such as the v' r v' g correlation, which arise 
from organized tilts within the convective structures, 
especially in the downflow plumes (e.g. Brummell et 
al. 1998, Miesch et al. 2000). 

In Figure 11 we show the components of J- r and J-g 
for cases A, AB, B and C , having integrated along 
co-latitude and radius respectively to deduce the net 
fluxes through shells at various radii and through 
cones at various latitudes, namely in the manner 



T r {r,6)r 2 sin Odd 



Mr, < 



i r sin 



>dr, 



(10) 



(11) 



1-bot 



and then identify in turn the contributions from vis- 
cous (V), Reynolds stresses (R) and meridional cir- 
culation (M) terms. This representation is helpful in 
considering the sense and amplitude of the transport 
of angular momentum within the convective shells by 
each component of J- r and J-g. 
Turning first to the radial fluxes in the leftmost of 



each pair of panels in Figure 11, we note that the in- 
tegrated viscous flux T r y is negative (where for sim- 
plicity we drop I), implying a radially inward trans- 
port of angular momentum. This property is in agree- 
ment with the positive radial gradient in the angular 
velocity profiles achieved in our four cases, as seen 
in Figure |I| in the radial cuts for different latitudes 
of 0. Such downward transport of angular momen- 
tum is well compensated by the two other terms J 7 r ,R 
and T r ,Mi having reached a statistical equilibrium of 
nearly no net radial flux, as can be seen by noting 
that the solid curve J- r is close to zero. Although all 
of our solutions possess complicated temporal vari- 
ations, our sampling in time to obtain the averaged 
fluxes suggest that we are sensing the equilibrated 
state reasonably well. As the level of turbulence is in- 
creased in going from case A to C, T r y reduces in am- 
plitude and the transport of angular momentum by 
the Reynolds stresses and by the meridional circula- 
tion change accordingly to maintain equilibrium. The 
meridional circulation as J r r ,M involves a strong domi- 
nantly outward transport of angular momentum. The 
Reynolds stresses as J~ r ,R vacillate in their sense with 
depth, though consistently possess outward transport 
in the upper portions of the domain. Case AB is dis- 



tinguished by J- r jt being directed outward through- 
out the domain. Detailed examination with radius 
and latitude of the Reynolds stress contributions to 
the angular momentum fluxes in equations (7-9) re- 
veals that the 'flux streamfunctions' (not shown) pos- 
sess multi-celled structures with radius at latitudes 
above 45° for all cases except AB. This striking dif- 
ference in case AB of having a big positive J 7 r ,R, ap- 
pears to influence the redistribution of angular mo- 
mentum at high latitudes. This may be key in the 
monotonic decrease of f2 with latitude of case AB ex- 
tending into the polar regions, and provides our first 
clue for how Issue 1 is resolved within this case. In 
a broader sense in considering all of our cases, we 
deduce that in the radial direction the transport of 
angular momentum is significantly affected by both 
the meridional circulation and the Reynolds stresses. 

The latitudinal transport of angular momentum 
J-g, in the rightmost of the panels in Figure 11, in- 
volves more complicated and sharper variations in 
latitude. This comes about due to the more intri- 
cate latitudinal structure of the differents terms con- 
tributing to the transport. Here the transport of an- 
gular momentum by Reynolds stresses J r g,R appears 
to be the dominant one, being consistently directed 
toward the equator (i.e. negative in the south hemi- 
sphere and positive in the north hemisphere). This 
is an important feature, since it implies that the 
equatorial acceleration observed in our simulations 
is mainly due to the transport of angular momentum 
by the Reynolds stresses, and thus is of dynamical 
origin. As we try to understand Issue 2, concerned 
with retaining a significant AO as the flow complex- 
ity is increased, we find that the variation of angular 
momentum fluxes by Reynolds stresses with increas- 
ing complexity along Paths 1 and 2 are fairly similar 
in character. Along both these paths the Reynolds 
stress fluxes remain prominent, and this appears to 
sustain the large AO, thereby resolving Issue 2 for so- 
lutions with the level of turbulence attained in cases 



C and D (the latter is not shown in Figure 11, but its 
transport properties are comparable to those of case 
C). Further, we see that the transport by merid- 
ional circulation Tg.u is opposite to Fg,R, with the 
meridional circulation seeking to slow down the equa- 
tor and speed up the poles. A distinguishing feature 
of case AB is that J-g m becomes small at latitudes 
above 30°, with the tendency of the meridional cir- 
culation to try to spin up the high latitudes sharply 
diminished compared to the other cases. This ap- 
pears to result from the strong meridional circulation 
in case AB being largely confined to the interval from 
the equator to 30° in latitude (Fig. ||a), with only 
a weak response at higher latitudes. This property 
of Tg m, together with the uniformly positive T r R, 
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Fig. 11. — Time average of the latitudinal line integral of the angular momentum flux T r (left panels in each pair) and of the 
radial line integral of the angular momentum flux Te (right panels) for cases A (top left), AB (top right), B (bottom left) and C 
(bottom right). The fluxes have been decomposed into their viscous (labelled V), Reynolds stress (R), and meridional circulation 
(M) components. The solid curves represent the total fluxes and serve to indicate the quality of stationarity achieved. The positive 
values represent a radial flux that is directed outward, and a latitudinal flux directed from north to south. The fluxes for A, AB, 
B and C have been averaged over periods in turn of 295, 275, 275 and 175 days. 



provides the second clue for how Issue 1 appears to 
be resolved by case AB. As the level of turbulence 
is increased, we find a reduction in the amplitudes 
of all the components of Fg, with T$y always being 
the smallest and transporting the angular momentum 
poleward in the same sense as Te,M- F° r ^e,R, this 
lessening amplitude appears to come about from the 
increasing complexity of the flows implying smaller 
correlations in the Reynolds stress terms, but it is 
likely that strengthening coherent turbulent plumes 
can serve to rebuild such correlations (Brummell et 
al. 1998). 

Our estimates of the latitudinal transports of an- 
gular momentum yield fairly good equilibration for 
cases A and AB, with little net latitudinal flux, but 
the more turbulent cases such as C are sufficiently 
complex that achieving such latitudinal balance is a 
slow process in the temporal averaging. We conclude 



that the Reynolds stresses have the dominant role in 
achieving the prograde equatorial rotation seen in our 
simulations, with its effectiveness limited by the op- 
posing transport of angular momentum by the merid- 
ional circulation. The viscous transports are becom- 
ing more negligible as we achieve more turbulent flows 
by reducing the eddy diffusivities. 

4.2. Baroclinicity and thermal winds 

Convection influenced by rotation can lead to lat- 
itudinal heat transport in addition to radial trans- 
port, thereby producing latitudinal gradients in tem- 
perature and entropy even if none were imposed by 
the boundary conditions. This further implies that 
surfaces of constant mean density and mean pres- 
sure will not coincide, thereby admitting baroclinic 
terms in the vorticity equations (Pedlosky 1987, Zahn 
1992). Baroclinicity has been argued to possibly have 







Fig. 12. — Temporal and longitudinal average for case AB of (a) the longitudinal velocity v<t,, (b) its derivative along the z 
axis, dv^/dz, (c) the baroclinic term in the meridional force balance (see equation 14), and (d) the difference between the last two 
terms [namely (b) minus (c)]. The results have been averaged over a period of 275 days. Panel (d) shows that there are major 
departures from a simple thermal- wind balance, especially near the top of the domain. The same color scale is used in panels (b), 
(c) and (d). 



a pivotal role in obtaining differential rotation profiles 
whose angular velocity, like the sun, are not constant 
on cylinders (e.g. Kitchatinov & Riidiger 1995). We 
shall here analyze our cases AB and C from that 
perspective, finding that though a small latitudinal 
entropy gradient is realized, the resulting differential 
rotation as exhibited in our solutions by the mean 
longitudinal velocity cannot be accounted for prin- 
cipally by the baroclinic term. To make such inter- 
pretation specific, we should turn as in Elliott et al. 
(2000) to the mean (averaged in longitude and time) 
zonal component of the curl of the momentum equa- 
tion (2), expressed as 



_d_ 



dvj_ 

dxi. 



d 
dxi 



V 



2Q 



dv<p 
dz 



(12) 



+ P 
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where the Einstein summation convention has been 
adopted, e represents the permutation tensor, and 
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is the derivative parallel to the rotation axis. This 
vorticity equation is helpful in examining the relative 
importance of different forces in meridional planes; 
here terms arising from Reynolds and viscous stresses 
are on the left and from Coriolis and baroclinic ef- 
fects on the right. If one were to simply neglect the 
Reynolds and viscous stresses, we obtain the simplest 



version of a 'thermal-wind balance' in which depar- 
tures of zonal winds from being constant on cylinders 
aligned with the rotation axis are accounted for by the 
baroclinic term involving crossed gradients of density 
and pressure, namely 



2 ^ 

dz 



'VpxVp 



With the superadiabatic gradient expressed as 

—vs = -4vp- \v p , 

c P IP p 



(13) 



(14) 



where 7 is the logarithmic derivative of pressure with 
respect to density at constant specific entropy, we can 
rewrite equation (|l3|) as 



dv<j> 1_ 

dz 2ft pcp 



VSxVp 



9 



dS_ 

2n o rc P ~d9 ' 



(15) 



having neglected turbulent pressure. Thus breaking 
the Taylor-Proudman constraint that requires rota- 
tion to be constant on cylinders, with dv^/dz zero, 
can be achieved by establishing a latitudinal entropy 
gradient. However, Reynolds and viscous stresses can 
also serve to break that constraint, and indeed we 
next show that those terms are at least as important 
as the baroclinic term. 



We turn in Figure 12 to an analysis of case AB in 
terms of how well is a simple thermal-wind balance 



achieved or violated. Figures 12 a, b display the tem- 
poral mean zonal velocity va, and its gradient dv^/dz, 
with the latter having pronounced variations at mid 
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Fig. 13. — As in Fig. [l^, analyzing the role of baroclinicity in the more turbulent case C in maintaining the differential 
rotation. There are significant departures from a thermal-wind balance in thin regions near the top and bottom of the shell. 



latitudes near the top of the spherical shell and oth- 
ers at lower latitudes near the bottom of the domain. 
The baroclinic term [as on right of equation (|l~5|)1 is 
shown in Figure |l^c, possessing the largest ampli- 
tudes close to the base of the shell at low latitudes, 
with a tongue connecting to mid latitudes in travers- 
ing the shell. The difference between this baroclinic 
term and the actual dv^/dz, as shown in Figure ^%(L 
is a measure of the effectiveness of a thermal-wind 
balance in case AB. It is evident that baroclinicity 
yields a fair semblance of a balance over much of the 
deeper layer, with the baroclinic term (Fig. |i~2|c) typ- 
ically being greater in amplitude than dv^/dz (Fig. 
|j~2]fr) there. However, the major regions of depar- 
ture with opposite signs in the two hemispheres show 
that in the upper domain, between latitudes of about 
15° and 45°, that balance is quite severely violated: 
there the Reynolds stress terms in equation ( |i~2| ) in- 
volving vortex tube stretching and tilting become the 
main players. This broad site coincides with regions 
of strong latitudinal gradient in v^, and is centered 
in latitude where the relative rotation changes sense 
from prograde to retrograde. What we have learned 
from this is that whereas the convection does estab- 
lish a latitudinal gradient of entropy that is needed 
for baroclinic terms to achieve aspects of thermal- 
wind balance over the deeper portions of the domain, 
the Reynolds stresses have an equally crucial role in 
the meridional force balance over portions of the up- 
per domain. The more turbulent case C is likewise 



analyzed in Figure 13, and it generally exhibits com- 
parable behavior. The baroclinic term (Fig. [Hfc) 
captures much of the dv^/dz variation (Fig. |l3|fr) at 
mid latitudes over most of the deep shell, but there 



are large departures (Fig. |l~3p) in thin domains near 
the top and bottom of the shell, again between 15° 
and 45° in latitude. Thus here too the Reynolds stress 
terms are significant players in the overall balance. 

The latitudinal entropy and temperature gradients 
established within our simulations should be exam- 
ined further. We show in Figure [L4| the time and lon- 
gitude averaged specific entropy fluctuations S and 
temperature fluctuations T for cases AB and C, pre- 
senting both color contour renderings across the shell 
and their variations with latitude at three depths. 
Our model AB, which exhibits the strongest differen- 
tial rotation, also possesses the greatest temperature 
and entropy contrasts with latitude. We see from 
the latitudinal cuts of temperature that a Af2 of or- 
der 30% involves a pole-equator temperature varia- 
tion of about 4 to 8 degree K, the pole being warmer. 
These temperature contrasts are very small compared 
to the mean temperature near the top of our domain 
of about 10 5 K, and of 2 x 10 6 K near its base. There 
is some evidence of a latitudinal variation in the pho- 
tospheric temperature of at least a few K with the 
same sense obtained from observations of the solar 
limb (e.g. Kuhn 1998), though relative variations of 
such small amplitude are very difficult to measure. 
We note that our temperature fields show some band- 
ing with latitude near the top of the domain, with the 
equator slightly warm, then attaining relatively cool 
values with minima at about latitude 35°, followed by 
rapid ascent to warm values at high latitudes. The 
behavior is monotonic with latitude at greater depths, 
as it is consistently so for entropy at all depths. These 
differences between temperature and entropy 
counted for by effects of the pressure field necessary 



21 




umtnflo (LtBjjj LitUudi fargl 

Fig. 14. — Temporal and longitudinal averages for cases AB and C of the specific entropy (upper panels) and temperature 
fluctuations (lower panels), accompanied by latitudinal profiles at the base {dash dotted line), at the middle {dashed) and at the 
top {solid) of the convective domain. The results have been averaged over periods in turn of 275 and 175 days. The presence of a 
latitudinal variation of entropy is consistent with the baroclinic term (shown in Figs. [l2| and |l^) and involves an equator to pole 
temperature contrast of at most 4 to 8 K near the top where the mean temperature is about 10 5 K. 



to drive the meridional circulation. 

In summary, although our solutions attain close to 
a thermal-wind balance over large portions of the do- 
main, the departures elsewhere are most significant. 
These arise from the Reynolds stresses that have a 
crucial role in establishing the differential rotation 
profiles realized in our simulations. The baroclinicity 
in our solutions, resulting from latitudinal heat trans- 
port that sets up a pole-to-equator temperature and 
entropy contrast, contributes to f2 not being constant 
on cylinders, but it is not the dominant player as envi- 
sioned in some discussions of mean-field models of so- 
lar differential rotation (e.g. Kitchatinov & Riidiger 
1995, Rekowski & Riidiger 1998; Durney 1999, 2000). 

5. CONCLUSIONS 



Our five simulations studying the coupling of tur- 
bulent convection and rotation within full spherical 
shells have revealed that strong differential rotation 
contrasts can be achieved for a range of parameter 
values. With these new models, we have focused on 
two fundamental issues raised in comparing the so- 
lar differential rotation deduced from helioseismology 
with the profiles achieved in the prior 3-D simulations 
of turbulent convection with the ASH code (Miesch 
et al. 2000, Elliott et al. 2000). As Issue 1, the sun 
appears to possess remarkably slow poles, with de- 
creasing steadily with latitude even at mid and high 
latitudes (Fig. |l|). In contrast, the previous models 
showed little variation in Q at the higher latitudes, 
having achieved most of their latitudinal angular ve- 
locity contrast Af2 in going from the equator to about 
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45°. As Issue 2, there was a tendency for AQ to di- 
minish or even decrease sharply within the prior sim- 
ulations as the convection became more turbulent, 
yielding values of A $7 that were becoming small com- 
pared to the helioseismic deductions. In seeking to 
resolve these two issues, we have explored two paths 
in parameter space that yield complex and turbulent 
states of convection. Path 1 involves decreasing the 
Prandtl number in the sequence of cases A, B and 
C, while keeping the Peclet number nearly constant. 
Path 2 maintains a constant Prandtl number as both 
the Reynolds and Peclet number are increased in the 
sequence of cases AB, B and D. On both paths the 
convective Rossby number has been chosen to be less 
than unity, thereby maintaining a strong rotational 
influence on the convection even as the flows become 
more intricate. 

In dealing with Issue 1, our case AB provides the 
first indications that it is possible to attain solutions 
in which the polar regions rotate significantly slower 
than the mid latitudes (Fig. |||) . There is a monotonic 
decrease from the fast (prograde) equatorial rate in 
to the slow (retrograde) rate of the polar regions. 
Further, that case AB has O nearly constant on radial 
lines at the higher latitudes, again in the spirit of the 
helioseismic inferences. We do not fully understand 
why in case AB such a strikingly different f2 profile re- 
sults, compared to that in our other solutions (and of 
the progenitor simulations) in which the contrast Aft 
is mainly achieved in the lower latitudes. Our prin- 
cipal clues come from Figure pi] where we find that 
only in case AB is the Reynolds stress component of 
the net radial angular momentum flux J- r R (through 
shells at various radii) uniformly directed outward. 
From having examined in detail angular momentum 
flux streamfunctions (not shown) with radius and lat- 
itude consistent with equations (7-9), we observed 
that the Reynolds stress contributions to such trans- 
port possessed multi-celled structures with radius at 
high latitudes in all the cases except AB. The single- 
cell behavior there for case AB appears to enable 
more effective extraction of angular momentum by 
Reynolds stresses from the high to the low latitudes, 
thereby yielding a distinctive rotational slowing of the 
high latitudes. Further, case AB possesses strong 
meridional circulations at low latitudes, but only fee- 
ble ones at latitudes above 30°, unlike other solutions 
such as case C (Figs. ||, ||). This yields a weak merid- 
ional component Fq.m (seeking to spin up the poles) 
to the latitudinal angular momentum flux, thereby 
allowing the equatorward transport by the Reynolds 
stress component Fq ; r to succeed in extracting angu- 
lar momentum from the higher latitudes. Such polar 
slowing also leads to case AB possessing the greatest 
AQ attained in our five simulations (Table 2). 



We also considered the possibility that the slow 
pole behavior in case AB may have baro clinic ori- 
gins. This can result from suitable correlations in 
velocity and thermal structures yielding a latitudi- 
nal heat flux which may produce substantial entropy 
variations at the higher latitudes, thereby leading to 
greater baroclinic contributions in equation (11) that 
defines the variation of mean zonal velocity. Exami- 
nation of Figure |l^ at high latitudes does not reveal 
a prominent baroclinic contribution, and this is con- 
sistent with the bland variation of entropy for case 
AB (Fig. U) at latitudes above about 40°. We con- 
clude that the origin of the slow rotation rate in po- 
lar regions appears to be primarily dynamical, be- 
ing associated with the Reynolds stress transports, 
and not with baroclinicity that arises from latitudinal 
heat transport serving to establish a sufficiently warm 
pole. Although case AB provides a solution that re- 
solves Issue 1 , it is unique in achieving this among our 
five simulations. It may be that in parameter space 
there only exists a small basin of attraction for such 
behavior, though we think it more likely that sev- 
eral solution states may coexist for the same control 
parameters, one of which exhibits the gradual rota- 
tional slowing at high latitudes, and others having 
most Q variations confined to low and mid latitudes. 
We plan to examine whether the slow pole character- 
istics of case AB can be maintained at nearby sites 
in parameter space if started from initial conditions 
corresponding to AB, and plan to report on this in 
the future. 

Issue 2 concerns sustaining a strong differential ro- 
tation with latitude as the convection becomes more 
complex. The two paths that we have explored in 
parameter space to achieve more complex and tur- 
bulent states yield relative angular velocity contrasts 
An/n in latitude that are comparable to values de- 
duced from helioseismology, with both our models 
and apparently the sun possessing a contrast of order 
30%. Further, this is accomplished while imposing an 
upper thermal boundary condition that ensures a uni- 
form emerging heat flux with latitude, as suggested 
in Elliott et al. (2000). Path 1 involving a decreasing 
Prandtl number is somewhat more effective in attain- 
ing large AQ as the solutions become turbulent than 
Path 2 which has the Prandtl number fixed at 0.25 as 
both diffusivities are decreased. This holds out hope 
that even more turbulent solutions will act likewise. 

We have shown that the strong AO, results from 
the role of the Reynolds stresses in redistributing the 
angular momentum. This transport is established 
by correlations in velocity components arising from 
convective structures that are tilted toward the rota- 
tion axis and depart from the local radial direction 
and away from the meridional plane. These yield 
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both v T v$ and vqv$ correlations necessary for the 
Reynolds stress contributions to the radial and lat- 
itudinal angular momentum fluxes analyzed in Fig- 
The fast downflow plumes have a dominant 



11 



ure 

role in such Reynolds stresses, much as seen in lo- 
cal studies (Brummell et al. 1998). Our simulations 
have attained a spatial resolution adequate to begin 
to attain coherent structures amidst the turbulence, 
which is believed to be a key in sustaining strong 
Reynolds stresses at higher turbulence levels. This 
has the consequence that all our spherical shells pos- 
sess fast prograde equatorial rotation relative to the 
reference frame. There are some contributions to- 
ward maintaining the differential rotation from the 
latitudinal heat transport inherent in our convection 
that serves to establish a warm pole (with a contrast 
of a few K) relative to the equator, with baroclinicity 
and a partial thermal-wind balance helping to yield 
equatorial acceleration. The meridional circulations 
generally work to oppose such tendencies by redis- 
tributing angular momentum so as to try to spin up 
the poles. Our simulations on Paths 1 and 2 confirm 
that strong differential rotation with fast equators has 
its primary origin in angular momentum transport as- 
sociated with the Reynolds stresses. Such prominent 
transports serve to resolve Issue 2. Our next chal- 
lenge is to satisfy Issue 1 simultaneously with Issue 2 
in the more turbulent solutions, which may also lead 
to f2 being more nearly constant on radial lines at 
mid to high latitudes. 

Although our results for f2 have made promising 
contacts with helioseismic deductions about the state 
of solar differential rotation in the bulk of the con- 
vection zone, there are also major issues that we 
have not yet tackled. We must evaluate more ad- 
vanced subgrid-scale terms in representing the unre- 
solved turbulence within such simulations, especially 



in the near-surface regions. Foremost are also ques- 
tions of how does the presence of a region of pene- 
tration below the convection zone influence the an- 
gular momentum redistribution in the primary zone 
above, and does the tachocline of shear that is es- 
tablished near the interface with the deeper radia- 
tive interior modify properties within the convection 
zone itself. We are keen to also investigate aspects 
of the rotational shear evident close to the solar sur- 
face. This is just now becoming computationally fea- 
sible, and involves extending our computational do- 
main upward and beginning to resolve supergranu- 
lar motions there, as discussed in DeRosa & Toomre 
(2001) in preliminary studies with thin shells. We are 
still at early stages with our simulations using ASH 
to study turbulent convection in spherical shells, yet 
it is comforting that the mean differential rotation 
profiles realized in some of the simulations are begin- 
ning to capture many of the dominant features for Q 
deduced from the helioseismic probing. 
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Pittsburgh Supercomputing Center (PSC). Much of 
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Table 1 

Parameters for the Five Simulations 



Case 


A 


AB 


B 


C 


D 




64, 128, 256 


64, 128, 256 


64, 256, 512 


192, 256, 512 


192, 512, 1024 


Ra 


3.1 xlO 4 


3.4 xlO 4 


1.4 xlO 5 


3.1 xlO 5 


6.5 xlO 5 


T a 


7.7 xlO 4 


3.1 xlO 5 


1.2 xlO 6 


5.4 xlO 6 


6.5 xlO 6 


Pr 


1 


0.25 


0.25 


0.125 


0.25 


Roc 


0.645 


0.662 


0.673 


0.682 


0.633 


V 


5.5 xlO 12 


2.8 xlO 12 


1.4 xlO 12 


6.8 xlO 11 


6.0 xlO 11 


K 


5.5 xlO 12 


1.1 xlO 13 


5.5 xlO 12 


5.5 xlO 12 


2.4 xlO 12 


Re 


28 


85 


170 


385 


410 


R 


0.10 


0.16 


0.15 


0.17 


0.16 


Pe 


28 


21 


43 


48 


103 



All simulations have an inner radius r^ot = 5.0 x 10 cm, an outer radius rt op = 6.72 x 10 cm, with L = 1.72 x 
10 10 cm the thickness of the computational domain. The number of radial, latitudinal and longitudinal mesh 
points are N r , Nq, N^. Here evaluated at mid-layer depth are the Rayleigh number R a = (—dp/dS)ASgL 3 /pun, 
the Taylor number T a = A&L^/v 2 , the Prandtl number P r = v/k, the convective Rossby number R oc = 
^/R a /T a P r , the rms Reynolds number R e = vL/u, the rms Peclet number P e = R e P r = vL/k, and the rms 
Rossby number R Q = uj/20, ~ v/2VLL, where v is a representative rms convective velocity. A Reynolds number 
based on the peak velocity at mid depth would be about a factor 4 larger. The eddy viscosity v and eddy 
conductivity k at mid depth are quoted in cm 2 s _1 . 



Table 2 

Representative Velocities, Energies and Differential Rotation 



Case 



Mid Convective Zone 



V<f, 



KE 



DRKE 



Volume Average 
CKE 



MCKE 



1.9xlO b (70%) 
2.3xl0 6 (36%) 
3.1xl0 6 (48%) 
4.3xl0 6 (54%) 
4.2xl0 6 (65%) 



1.0 x 10 4 (0.37%) 
(0.32%) 
(0.38%) 
(0.42%) 
(0.46%) 



A 

AB 
B 
C 
D 



46 
50 
57 
68 
72 



40 
47 
56 
67 
67 



69 

124 

115 

122 

108 



44 
53 
59 
70 
64 



92 

142 

140 

155 

146 



74 

87 

99 

117 

111 



2.7x10° 
6.5xl0 6 
6.5x10° 
7.9xl0 6 
6.5xl0 6 



8.2x10° 
4.2xl0 6 
3.4x10° 
3.6x10° 
2.3x10° 



(30%) 
(64%) 
(52%) 
(46%) 
(35%) 



2.1x10" 
2.5xl0 4 
3.3xl0 4 
3.0xl0 4 



In the five cases, temporal averages at middayer depth in convection zone of rms components of velocity v r , vg, 
Vfj, and of speed v, and of fluctuating velocities v'^ and v' (with temporal and azimuthal mean subtracted), all 
expressed in ms" 1 . Also listed are the time averages over the full domain of the total kinetic energy, KE, that 
associated with the (axisymmetric) differential rotation, DRKE, with the (axisymmetric) meridional circulation, 
MRKE, and with the non-axisymmetric convection itself, CKE, all in ergcm~ 3 . The relative latitudinal contrast 
of angular velocity AQ/Q between 0° and 60° near the top of the domain are stated. 



